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SUMMARY 


/ a 3 ^ ^ 

A method of analyzing flow through a turbomachine is presented that is 
suitable for computer programing. It is assumed that a mean stream surface 
from hub to shroud between blades is known. On this stream surface a two- 
dimensional solution for the velocity and pressure distributions is obtained, 
and then an approximate calculation of the blade surface velocities is made. 
This method is based on an equation for the velocity gradient along- an arbi- 
trary quasi -orthogonal rather than the normal to the streamline as used in pre- 
vious methods. With this new method a solution can be obtained in a single 
computer run, even for cases where the distance between hub and shroud is great 
and there is a change in direction from radial to axial within the rotor. The 
method was successfully applied to a turbine with this type of geometry. These 
results are given as a numerical example, and the Fortran computer program is 
included. 



INTRODUCTION 

Quasi-three-dimensional methods have been developed for analyzing flow 
through mixed-flow turbomachines. One such method (ref. l) is based on the as- 
sumption of axial symmetry and on an equation for the velocity gradient along 
the normal to the projection of the streamlines on a plane containing the axis 
of rotation. This basic method was used (ref. 2) to redesign the hub-shroud 
profile of a compressor rotor. The results of reference 3 showed improved per- 
formance for impellers redesigned by this method. A computer program using 
this method for the design of pump impellers is described in reference 4. In 
reference 4, the same velocity gradient equation given in reference 1 is de- 
veloped without the assumption of axial symmetry but with the assumption of a 
known stream surface that extends from hub to shroud. The examples used in the 
aforementioned references were all compressors or pumps, but the method is 
equally applicable to turbines. 

These methods use streamlines and their normals to establish a grid for 



the solution. In cases where the distance between hub and shroud is great and 
there is a large change in flow direction within the rotor, however, the nor- 
mals vary considerably in length and direction during the course of the calcu- 
lations. Therefore, it becomes difficult to obtain a direct solution on the 
computer without resorting to intermediate graphical steps. The use of nor- 
mals, however, is not essential to the method, and it appeared possible to 
overcome this difficulty by the use of a set of arbitrary curves from hub to 
shroud instead of streamline normals. These arbitrary curves will be herein- 
after termed quasi -orthogonals. The quasi -orthogonals are not actually ortho- 
gonal to each streamline, but merely intersect every streamline across the 
width of the passage. The quasi -orthogonals remain fixed regardless of any 
change of streamlines. By using this technique, it appeared possible to de- 
velop a computer program that would calculate the velocity and pressure distri- 
butions without any intermediate graphical procedures even for turbomachines 
with wide passages and a change of direction from radial to axial within the 
rotor blade. 

In view of these considerations, a method of analysis utilizing quasi- 
orthogonals in lieu of streamline normals was developed. This report presents 
the analysis method and contains a discussion of the numerical techniques re- 
quired for obtaining solutions with a digital computer. The computer program 
developed during this study is included. As a numerical example of the appli- 
cation of the analysis method, a radial-inlet mixed-flow gas turbine of high 
specific speed is analyzed. Such a turbine, which may have application in gas 
turbine cycle space power systems, has a rotor-channel geometry for which this 
method, as compared to previous methods, can yield a quick and direct solution. 


METHOD OF ANALYSIS 

The analysis to be presented herein is basically the same as those pre- 
sented in references 1, 2, and 4. As pointed out In the INTRODUCTION, the 
major difference is the use of fixed arbitrary quasi -orthogonals rather than 
streamline normals to establish a grid for the solution. Another difference is 
that the reference analyses were based completely on the assumption of isen- 
tropic flow, while in this analysis a correction for a loss in relative total 
pressure is included in the continuity equation to account for blade losses. 

This analysis, as that of reference 4, is based on the assumption of a 
mean flow surface between blades. In general, this surface is assumed to be 
parallel to the mean blade surface, with arbitrary or empirical corrections 
made to take care of the difference between the flow angle and the blade angle 
at the inlet and at the outlet. One factor that is not accounted for by this 
assumption is that the actual mean stream surface twists considerably in a 
mixed-flow turbine. Despite this assumption, however, reference 3 shows the 
value of the analysis method by the improved performance of compressors re- 
designed in accordance with this assumption. Reference 5 shows that a two- 
dimensional solution for a particular compressor, when compared to a three- 
dimensional solution, gives values of the through-flow component of velocity 
that are of sufficient accuracy for engineering analysis. For convenience, the 
mean stream surface is projected on a plane containing the axis of rotation. 
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This plane is called the meridional plane. The projections of the streamlines 
on the meridional plane are called meridional streamlines. 


Analytical Equations 


Equations (l) and (2) give the velocity gradient along an arbitrary quasi - 
orthogonal in the meridional plane 


dW 

ds 


= (a + B + C ^ + D ¥ + fe - 0) 

\ ds ds/ ds ds yds ds f 


CD 


where 


. cos a cos^B sin^R , . 

A = — - + sm a sm 3 cos 3 


be 


^ sin a cos 2 3 , . „ bG 

B = - + sm a sm 3 cos 3 V” 

r n oz 




dW, 


C = sin a cos (3 — — - 2<jd sin 8 + r cos 0l — — + 2od sin a 
dm 


VfWe 

\dm 


dW 


D = cos a cos 3 — — + r cos 3 
dm 


/dWg 

\dm 


ydm 

+ 2a> sin 


) 




in 


de 

SI 


( 2 ) 


The coordinate system and the notation are shown in figures 1 and 2. (All 
symbols are listed in appendix A.) Equation (l) is derived in appendix B. 



Figure 1. - Coordinate system and velocity components. 


The value of the parameters hj_ and 
A associated with a point inside the rotor 
is the value of that parameter at the inlet 
for the streamline which passes through 
that point. Then dh^/ds refers to the 
total enthalpy at the inlet as a function 
of the distance along the arbitrary meri- 
dional quasi -orthogonal near the point in 
question. 


X \-Streamline 



Figure 2. - Component of relative velocity W n normal to 
arbitrary quasi-orthogonal. 
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Radius from axis of rotation, 



Axial coordinate, z, in. 

Figure 3. - Profile of rotor in numerical example. 


In this analysis , the 
arbitrary quasi -orthogonals 
were chosen to be straight 
lines from hub to shroud 
(see fig. 3). At the inlet 
and outlet, the lines were 
chosen as the leading and 
trailing edges, respectively. 


In addition to equa- 
tion (l), which is a force 
equilibrium equation, the 
continuity equation must be 
satisfied. This is done by 
requiring that the calculated 
weight flow across any line 
from hub to shroud be equal 
to the specified weight flow 
through the turbomachine . 

For this the density must be 
known. If the velocity is 
known, the density may be 
calculated by equations (3) 
to (5) following. Equa- 
tion (B9) is 
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which gives the static density at any point once the velocity is known if inlet 
total conditions are specified. 


To account for losses, it is necessary to make a correction to the above 
calculated density. One way to do this is to assume a loss in relative total 
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( 5 ) 

This gives the static density with a specified loss in relative total pressure. 
The temperature ratios can be calculated from equations ( 3 ) and (4) . It is 
assumed that inlet total conditions are known. 

Weight flow across a quasi -orthogonal can now be computed by 

w = N f pW n r A0 ds (6) 

•'O 

where A 6 is the angular distance between blades and W n is the component of 
W normal to the surface of revolution generated by the fixed line. Prom 
figure 2, it can be seen that 


W n = W m cos 0+ “ a ) 


To get A 6, use is made of the fact that 
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where tg is the tangential thickness. If the thickness normal to the mean 
blade surface t n is specified. 
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Note that here d0/dr and de/d z refer to the mean blade shape and not the 
assumed mean stream surface. Equations (3) to (5) and (7) to (9) give the nu- 
merical data for equation (6), which can be integrated by use of a spline fit 
approximation (see appendix C) . 

With the velocities on the mean stream surface calculated, blade surface 
velocities can be calculated by any of several approximate methods. One method 
that gives good results, when compared with a relaxation solution of the po- 
tential flow equation for a surface of revolution, is based on absolute irrota- 
tional flow and linear velocity distribution between blades. The following 
equations based on these assumptions are equations (16) and (17) of refer- 
ence 6. 


w t = 
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( 10 ) 


W z = 2W - W t 


The derivative can be evaluated by use of a spline fit curve. 

Equations (3) to (5) and the equation of state can be used to calculate 
the static temperature, density, and pressure on the blade surfaces. 


Numerical Techniques and Procedure 


The first step in the analysis is the numerical evaluation of the param- 
eters a, p, r c , bo/br, b6/bz, dr/ds, dz/ds, dW^dm, and dW^/dm for use in 
equations (l) and (2). In order to evaluate the parameters a, p, and r c a 
streamline geometry must be established. First fixed straight lines (quasi- 
orthogonals) are drawn from hub to shroud along which the velocity gradient for 
an assumed streajn surface will be determined. For an initial approximation to 
the streamlines, each quasi -orthogonal can be divided into a number of equal 
spaces, as shown in figure 3. The success of the method is based on the fact 
that, for a reasonable assumed streamline pattern, the geometrical streamline 
parameters involved are not too different from those of the final solution. 

Bv means of a spline fit approximation (see appendix C) , dr/dz and 
d^r/dz* can be determined at each of the points established. Then 
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= tan-1 


dr 

dz 


and 
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The reciprocal of the radius of curvature (the curvature) is computed to avoid 
division by zero in case d^r/dz^ = 0. 


For the remaining parameters, the mean stream surface 0 = 0( r,z) between 
blades is needed; it must be given in such a manner that bo/br and bo/bz 
can be determined at any given point. The spline fit curve can assist in this. 
When do/dr and be/bz are known, p can be calculated from 


tan P 


d0 /be dr be dz\ (be be > 

r d£ = r {^ to + 51 dm/ = r \3F Sln a + S? cos a 


( 12 ) 


For the initial calculation, W may be assumed constant throughout the rotor. 
From figure 1, it is seen that 


and 


W m - W cos p 


Wg = W sin p 

Since the distance along the meridional streamline m is known, dW^/dm and 
dWg/dm can then be determined by the spline fit curve. Since dr /as and 
dz/ds are determined by the angles of the quasi-orthogonals, all the quanti- 
ties necessary for the calculation dW/ds from equation (l), except W it- 
self, are now determined. 


The next step is the numerical integration of equation (l), which is in 
the form 


dW 

ds 


f(W,s) 


where f is known only for a finite number of values of s. For a given ini- 
tial velocity on, say the hub, the velocity distribution along the quasi- 
orthogonal can be approximated by 

vi - + (f \ 

where the subscripts denote the number of the streamline, and As is the dis- 
tance along the quasi -orthogonal between streamlines. For an improved esti- 
mate, a Runge-Kutta method can be used. The following is a particular Runge- 
Kutta method that is well adapted for this case. Let 
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This avoids an obvious bias due to using the derivative at the beginning of the 
interval (see fig. 4) and gives a higher order approximation. For a mathema- 
tical analysis and error estimate, see reference 7. For the calculation of the 
quantity (dW/ds)j, equation (l) is used with the parameters calculated for the 
j th streamline and Wj . To calculate (dW/ds)j+i, the parameters calculated for 
the (j + l) st streamline are used and Wj+q is used for the velocity W in 
equation (l) . 

It should be noted that this method of integrating equation (l) involves 
much less computation than solving equation (l) directly and then numerically 
evaluating the resulting integral (e.g. , eq. (9) in ref. l) . This is espe- 
cially helpful for hand computation and is also helpful in simplifying computer 
programing. Accuracy is probably comparable; the method used here certainly 
gives satisfactory accuracy if the streamlines are spaced closely enough so 
that the velocity does not vary more than about 30 percent between streamlines. 
In the numerical example, the results using five streamlines did not differ 
appreciably from those using twenty streamlines. 


Completing this computation for a quasi -orthogonal from hub to shroud re- 
sults in the complete velocity distribution along that line based on the ini- 
tial estimate of the velocity on the hub. Equations (3) to (5) and (7) to (9) 
can be used to compute the integrand in equation (6). The numerical integra- 
tion can be performed by use of a spline fit approximation (see appendix C) . 
The computed total weight flow is then compared with the actual weight flow. 



Figure 4. - Approximation to solution of differential equation dW/ds - f(W, s). 


If the computed weight flow is too 
small, the velocity on the hub is 
increased, and vice versa. Then 
the velocity distribution and the 
weight flow are recalculated. The 
computed weight flow is a function 
of the assumed hub velocity; 
therefore, after two values of 
weight flow are computed, linear 
interpolation or extrapolation can 
be used to get an improved esti- 
mate for the hub velocity. A few 
iterations will determine the hub 
velocity that will give the cor- 
rect weight flow. 


— O — Obtained by Inverse 
interpolation 

— □ — Computed from eq. (7) 



Figure 5. - Weight flow distribution along quasi - 
orthogonal. 



for this together with the listing of 
ample are given in appendix D. 


From equation (6) the weight flow 
distribution along the quasi -orthogonal 
from hub to shroud can also be obtained. 
Inverse interpolation (by a spline fit 
approximation), can be used to determine 
the spacing of the streamlines on the 
quasi -orthogonal that will give equal 
weight flow between any two adjacent 
streamlines (see fig. 5). When this is 
done for every quasi -orthogonal from in- 
let to outlet, a new estimate for the 
meridional streamline pattern is ob- 
tained. This pattern, together with the 
calculated velocity distribution, can 
then be used for further iterations; 
however, using this estimate generally 
results in overcorrection. Therefore, 
only a fraction of the calculated cor- 
rection was made. Another problem is 
the tendency for the newly computed 
streamline to be less smooth than the 
previous streamline. If a computation 
is based on a set of streamlines that 
are not extremely smooth, the calculated 
streamline corrections become erratic. 
Thus it is important to be sure that the 
streamline estimate to be used for the 
following iteration be as smooth as pos- 
sible. Several methods of accomplishing 
this have been tried. The method that 
was successful for the cases tried was 
to use a streamline correction at each 
point of one-tenth the calculated cor- 
rection. With this the streamlines re- 
mained smooth, and a solution was 
reached in a single computer run, re- 
quiring about 50 iterations. Computer 
execution time was 2 minutes (on the 
IBM 7094). The computer program used 
omputed results for the numerical ex- 


NUMERICAL EXAMPLE 

The procedure outlined herein has been programed for solution on a digital 
computer. The following results were obtained for a particular turbine. The 
hub-shroud profile and quasi -or thogonals are shown in figure 3, together with 
the equally spaced streamlines used for the initial assumption. The blade has 
radial blade elements with the blade shape indicated in figure 6. There are 
13 blades, with no splitter blades. The rotational speed was 51,500 rpm, and 
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the fluid was air. The weight flow was 0.984 pound per second, inlet total 
temperature was 592° R, Vg at the inlet was 1010 feet per second, and the 
total inlet pressure was 42.5 pounds per square inch. The normal blade thick- 
ness was given by means of tabulated values on a grid. Blade thickness at any 
given point was obtained by linear interpolation. It was assumed that hj and 
A are both constant from hub to shroud. 

At the inlet, the flow surface was assumed to deviate from the blade sur- 
face in order to agree with the flow direction coming into the rotor. This 
angle at the inlet was -35°. The meridional streamlines are approximately ra- 
dial at the inlet, so that the stream surface was assumed to be independent of 
z where it deviates from the blade surface. The 6 coordinate was assumed to 
vary as the cube of r (and independent of z) for a given distance from the 
inlet. Let r^ denote the radius where the mean stream surface is assumed to 
deviate from the mean blade shape. Equation (13) of reference 6 gives an ap- 
proximate equation for determining r b , which may be written as follows: 


,-0.71 A0 


The equation of the stream surface for r > r b is 


(r - r b ) 3 tan 
3r i( r i “ r b) 2 


which, when differentiated, becomes 


de 

or 


(r - r b ) 2 tan 
r i( r i - r b) 2 


( 14 ) 




Figure 7. - Meridional projection of mean stream Figure 8. - Static pressure contours on mean stream Figure 9. - Relative velocity contours of mean stream 
surface for numerical example. surface for numerical example. surface for numerical example. 
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This is used in equations (2) and (12) 
when r > r b hut not in equation (9), 
since equation (9) refers to the blade 
shape. For the numerical example, r^ 
is about 1.60 inches. At the outlet it 
is assumed that the mean stream surface 
would follow the blade. There was also 
assumed to be a 2 . 5-pound-per-square- 
inch loss of total relative pressure, 
varying linearly from inlet to outlet. 

The calculated streamlines are 
shown in figure 7. Since the solution 
is restricted to the rotor blade, the 
streamlines near the outlet do not show 
the effect of downstream geometry. Some 
of the other calculated information is 
shown in the figures following. Fig- 
ures 8 and 9 show lines of constant 
pressure and constant relative velocity, 
while figure 10 shows blade loading dia- 
grams at the hub, the mean surface of 
revolution, and the shroud. 

Figure 8 shows that the pressure 
level is always decreasing in the direc- 
tion of flow. This is, of course, nor- 
mal for a radial turbine. In figure 9, 
it is seen that velocities are generally 
increasing, except along the hub near 
the inlet where they decrease slightly. 
Though not desirable, this can be toler- 
0 .2 .4 .6 .8 1.0 ated because of the favorable pressure 

Percent of distance along meridional streamline gradient. More serious are the negative 

(c) Shroud. velocities in the blade loading diagrams 

Figure 10. - Blade loading diagram for numerical example. (fig* 10) . This indicates an eddy on 

the trailing surface of the blade near 

the inlet , which would result in turbulence and mixing losses. Also^ a severe 
decreasing velocity gradient is indicated on the suction surface near the hub 
and near the shroud. This could lead to flow separation with accompanying high 
losses. 


CONCLUDING REMARKS 

A method of analysis of turbomachines is presented that is suitable for 
computer programing. The method and the results are similar to that obtained 
by other streamline analysis methods (e.g., refs. 1, 2 , and 4). The difference 
here is that velocity gradients are given along arbitrary quasi -or thogonals, 
rather than the normal to the streamlines as has been done in previous methods. 
The value of the method lies in the fact that a solution can be obtained in a 
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single computer run even for cases where the distance between hub and shroud 
is great and there is a change of direction from radial to axial within the 
rotor. The method was successfully applied to a turbine with this type of 
geometry. These results are given as a numerical example, and the Fortran com- 
puter program is included in appendix D. 

A more accurate hub-to-shroud analysis could be made by using information 
from a blade-to-blade streamline analysis. A blade-to-blade analysis would 
give a better approximation to the mean stream surface and also would give the 
blade-to-blade streamline spacing. Continuity would then be checked between 
the two hub-to-shroud stream surfaces instead of between blades. 


Lewis Research Center 

National Aeronautics and Space Administration 
Cleveland, Ohio, September 15, 1964 
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APPENDIX A 


SYMBOLS 

A parameter , eq. (2) 

a parameter, eq. (B14) 

B parameter, eq* (2) 

b parameter, eq. (B14) 

C parameter, eq. (2) 

c parameter, eq. (B14) 

stagnation speed of sound at inlet, ft/sec 
c p specific heat at constant pressure, (ft) (lb) /(slug) (°R) 

D parameter, eq. (2) 

f any function 

g acceleration due to gravity, ft/sec^ 

h static enthalpy, (ft) (lb) /slug 

m distance along meridional streamline, ft 
N number of blades 

n distance along normal to meridional streamline, ft 
p absolute static pressure, lb/sq ft 

Ap" loss in relative total pressure between inlet and any point 
q distance along arbitrary three dimensional curve, ft 

R gas constant, (ft) (lb) /(slug) (°R) 

r radius from axis of rotation, ft 

r b radius at which assumed stream surface is tangent to mean blade shape 
r c radius of curvature of meridional streamline, ft 
s distance along arbitrary quasi -orthogonal in meridional plane, ft 
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T temperature, °R 

t time, sec 

t n blade thickness normal to blade mean surface, ft 

t q blade thickness in circumferential direction, ft 

u unit vector 

V absolute fluid velocity, ft/sec 

W relative fluid velocity, ft/sec 

v weight flow crossing surface of revolution generated by quasi -orthogonal 
between hub and given point on quasi -orthogonal 

x x-coordinate 

y y-coordinate 

z axial coordinate 

a angle between meridional streamline and z-axis, radians 

(3 angle between relative velocity vector and meridional plane, radians 

T ratio of specific heat 

6 relative angular coordinate, radians 

A0 angle between blade surfaces at given point, radians 

A prerotation r^V^, sq ft/sec 

p mass density, slugs/cu ft 

cp absolute angular coordinate, radians 

\|/ angle between quasi -orthogonal and radial direction, radians 

ad rotational speed, radians/sec 

Subscripts : 

i inlet 

isen isentropic 

i number of streamline 
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I 


leading surface 


m component in direction of meridional streamline 

n normal component 

r radial component 

s shroud 

t trailing surface 

x x-component 

y y- component 

z axial component 

0 tangential component 

Superscripts : 

— vector quantity 

r absolute stagnation condition 

relative stagnation condition 



APPENDIX B 


DERIVATION OF THE VELOCITY GRADIENT EQUATION 

Euler's force equation for a nonviscous fluid is 

dV _ 1 _ 

dt p 


(Bl) 


This is simply expressed as three scalar equations in fixed rectangular coordi- 
nates x, y, and z. To reduce the problem to a steady-state condition, equa- 
tion (Bl) should be expressed in__ 
terms of the relative velocity W 
and the pressure gradient relative 
to a rotating cylindrical coordi- 
nate system r, 0, and z. The no- 
tation u* is used to denote a 
unit vector in the x direction; 
similar notation is used for the 
other coordinates. It should be 
y noted that the directions of the 
vectors u r and are functions 

of t as well as of cp. It is 
seen from figure 11 that 



(sin tf)u 


Figure 11. - Relations between unit basis vectors in absolute rectangular 
coordinates and relative cylindrical coordinates. 


u r - (cos cp)u x + (sin cp)uy 

U0 = - (sin + ( cos 9)uy 

Differentiating equation (B2) results in 
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Since 


VpU r + + ^z u z^ equation (B3) can be used to get 
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The pressure gradient will now be expressed in the relative coordinates. 
For the fixed cylindrical coordinates r, cp, and z. 

Note that actually u^ = ug, when cp and 6 refer to the same point (u~ 
varies with time, since q> varies with time for constant 0) . Also * 
dp/ckp = dp/d6, since bcp/ho = 1. This gives 

7p= ^^ + rSe {l0 + ll iI z (B5) 

Noting that W r = V r , W z = V 2 and Vq = \Jq + cor, substituting equa- 
tions (B4r) and (B5) in equation (Bl), and equating coefficients of u q, and 

u z result in 


dW r (W 0 + cnr) 2 
dt r 


1 

p or 


(B6a) 


1 d(rW g + or 2 ) ^ ^ 
r dt pr <5© 


(B6b) 


^z m _ 1 

dt p oz 


(B6c) 


Now an expression for the directional derivative of the relative velocity 
in any direction will be derived. The parameters in this expression require 
the knowledge of the streamline passing through a given point; however, once 
the streamline is known, the velocity gradient in any direction can be com- 
puted . 


If q denotes the distance along an arbitrary curve, the directional de- 
rivative of the pressure p along this curve is 

ci£ d£ dp dz 
dq “ 5r dq + ^0 dq + 5z dq 

Using equations (B6) gives 


1 d£ 
P dq 


dW r 

dt" 


(W 0 + cnr) 2 
r 



d(rW e + or 2 ) de 
dt dq 


dW 2 

dt 


dz 

dq 


(B7) 


Equation (B7) is an expression for the pressure gradient in the q direction. 
It is necessary to find a relation between the velocity gradient and the pres- 
sure gradient. This is easily done under the assumption that the flow is isen- 
tropic, so that 
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& = dh 
p 

Wow multiplying equation (B6a) by W r = dr/dt, equation (B6b) by W g = r d©/dt, 
and equation (B6c) by W z = dz/dt, then adding and combining terms yield 

1 dW 2 _ g dr 1 d£ _ d(r 2 ) dh 

2 dt " “ r dt ' p dt 2 dt ~ dt 

which is the energy equation for isentropic flow. Integrating from the inlet 
along a streamline results in 

W 2 - W 2 = 0 ) 2 (r 2 - r 2 ) - 2(h - h ± ) (B8) 

Since V m = W m and Vg = Wg + cur, 

V 2 _ a W 2 - w| = W 2 - Vg + 2^osr - oo 2 r 2 


or 


V 2 » W 2 + 2Vgcnr - cn 2 r 2 


hence , at the inlet, 


V? W? + 2cuA - 0 D 2 r? 

hp = hi + — = hi + g 

Substituting this for h i in equation (B8) gives 


h = h^ - cuA + 


cu 2 r 2 - W 2 


(B9) 


Since the flow is assumed isentropic, differentiating results in 

1 d£ _ dh dh i dA ? dr dW 
p dq ” dq ** dq ^ dq r dq dq 

Substituting this equation in equation (B7) yields 


dW 

dq 


1 ^i 


CD 


dA 


W dq W dq 


1 _ 

~ + W dt 


(Wg + cur) 2 


rW 


dr 1 d ( rW 6 + d© 


dq + W 


dt 


i. dz 

dq + W dt dq 


(BIO) 


Note that W m = W cos P and 
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da da dm __ w m 

dt dm dt r c 

Using this and differentiating W = W m sin a and W 7 = W m cos a result in 


__ W^cosVcos 


dW z cos^p si] 


+ W sin a cos p 


+ W cos a cos p 


(Bll) 


Also, 

1 d.(rV e ) ± d(rW 0 + r 2 ^) dW 0 

y — ^ ^ = r cos 3 + W sin a cos 0 sin 0 + 2rco sin a cos 0 


(B12) 


Using equations (Bll) and (B12) and the fact that Vq =* W 0 + oar, and 
Wq = W sin 0 in equation (BIO) gives 


— = a — + b — + c — + — ( _ a) — 

dq dq dq dq W \dq dq 


(B13) 


where 


W cos^p cos a W sin^p 


+ sin a cos p ;r— - 2 cjd sin p 


W cos p sin a 


+ cos a cos p 


(B14) 


C - w sin a cos p sin p + r cos p ^^ — + 2cjd sin 

The meridional plane analysis is concerned with the projection of the 
curve q onto the meridional plane. This projected curve will be the quasi- 
orthogonal. Letting s denote the distance along this meridional projection, 
then 


® = , h dz d£ 1 ( dh i dA 


ds dq * ds a ds + b ds + C ds + W i ds 


(B15) 


If the line s is a normal to the meridional streamline, then s = n, 
dr/ds = dr/dn = cos a, and dz/ds = dz/dn = - sin a, and equation (B15) reduce 
to equation (B24) of reference 4. 
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The quantities dr/ds and dz/ds in equation (B15) are determined by the 
parametric equations for the arbitrary curve in the meridional plane , r = r(s) , 
and z = z(s). The quantity d0/ds refers to the change in 0 in the actual 
curve q. If the curve q lies on a hub-to-shroud surface, which can be de- 
fined by 0 =* 0(r,z), then 


d£ _ dr + S0 dz 
ds dr ds 5z* ds 


(B16) 


By substituting equations (B14) and (B16), equation (BIS) can be rewritten in 
the following form: 


dW 

ds 


A a£ + B^W + C^+D|i 

ds ds / ds ds 



1 

W 


(i) 



Equation (l) is written in a form that is convenient for numerical solu- 
tion. 
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APPENDIX C 


USE OF SPLINE FIT CURVES 

If a set of function values corresponding to a set of arguments is given, 
there are several ways a curve can be fitted through these values so as to ap- 
proximate the original function with these values. The classical way is by an 
n th -degree polynomial for n + 1 points. This may not be satisfactory, how- 
ever, for a large number of points, especially for computing derivatives or 
curvature at end points. Another technique is to use fewer points to determine 
some sort of piecewise polynomial, but this does not lead to a smooth curve. A 
method that has received much attention recently is the piecewise cubic, with 
matching first and second derivatives, commonly referred to as a spline fit 
curve. Since for small slopes, the second derivative approximates the curva- 
ture of a function, the strain energy of a spline can be approximately mini- 
mized by minimizing /[f"(x)] 2 dx, where f(x) denotes the curve described by 
the spline. The spline fit curve has this property. This is proven in refer- 
ence 8. Thus the spline fit curve is a mathematical expression for the shape 
taken by an idealized spline passing through the given points. In reference 8, 
a simple procedure is outlined for determining the spline fit curve when the 
coordinates of the points are given together with two arbitrary end conditions. 
The end condition actually used in the computer program was that the second 
derivative at an end point is one-half the second derivative at the next point. 
This is equivalent to bending the spline beyond the last point slightly, in- 
stead of just letting it be straight. The spline fit curve provided a simple 
analytical method of determining many of the parameters in the equations. The 
spline fit curve was used to determine first and second derivatives, curvature, 
interpolated function values, interpolated derivatives, and for integration. 

One further point concerning the spline fit should be mentioned; that is, 
the approximation to an actual spline curve is dependent on the slope not being 
too large. Experimentally, good results are obtained if the absolute value of 
the slope is not greater than one. In applying this method to streamlines on a 
radial turbine, there is a problem since the angle may be around -90° at the 
inlet. This is easily overcome by rotating the coordinate axes 45° so that the 
maximum slope is about one. 
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APPENDIX D 


FORTRAN PROGRAM USED FOR NUMERICAL EXAMPLE 
Description of Main Program 

The FORTRAN program listed herein is the one used in the numerical ex- 
ample. It is written in FORTRAN IV and was run on an IBM 7094 digital com- 
puter. The program closely follows the steps given in the section on numerical 
procedure. The list of program variables preceding the program indicates the 
equation that is used to calculate a variable or the equation in which it is 
used. In the program, the number of the streamline is denoted by K and the 
number of the quasi-orthogonal by I. The inlet or the hub is denoted by 1. 

The program is written so that all linear measurements are in inches, angles are 
in degrees, and pressure is in pounds per square inch for both input and output. 
Units are changed to feet and radians for computation in the program. All other 
quantities are in the units specified in appendix A. 

It will be noted that a complete listing of input data cards is printed 
out. In the sample program, for example, the listing gives all the data used 
as input for the program. All input statements precede the comment card 
END OF INPUT STATEMENTS. 


Program Variables and Definitions 
A temporary storage 

AB(j) temporary storage 

AC(j) temporary storage 

AD(j) temporary storage 

AL(l,K) a 

ALM A (input variable) 

AR R (input variable) 

B temporary storage 

BA(K) total weight flow between hub and K th streamline 

BCDP integer (input variable); 1 will give DN, WA, Z, and R as output 

on cards in binary form after final iteration, for use as input 
for alternate conditions; 0 will cause this to be omitted 

BETA(I,K) p 
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BETAD 


BETAT 
BET IN 
C 

CAL(I,K) 

CBETA(I,K) 

CEF 

Cl 

CORFAC 

COSBD 

COSBT 

CP 

CURV(l,K) 

DELBTA(I) 

DELTA 

DEWSTY 

DN(l,K) 

DEDM( I) 

DTDR(I) 

DTDZ(I) 

DWMDM(I) 

DWTDM(I) 

E 

ERROR 

ERROR1 


$l> eq. (10) 

P t , eq. (10) 

P^ (input variable), eq. (14) 
temporary storage 
cos a 
cos (3 

tan p/r i (r j _ - r b ) 2 , eq. (14) 
c i 

percentage of calculated streamline correction to be used for next 
iteration (input variable) 

cos @ 2 > eq. (10) 

cos P b , eq. (10) 

C P 

l/r Q 

P t - Pj, eq. (10) 

calculated streamline correction (fig. 5) 

Pg 

distance along quasi -orthogonal from hub 

^ [(rto + W sin p)r A6] , eq. (lO) 

c)0/dr, eq. (2) and (l2) 

bd/bz, eq. (2), (9), and (l2) 

dW m / dm, eq. (2) 

dW^/dm, eq. (2) 

temporary storage 

maximum calculated streamline correction for present iteration 
(fig. 5) 

ERROR from previous iteration 
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EXPON 


G 

GAM 

HR 

HZ 

I 

USD 

ITER 


J 

K 

KMX 

KMXM1 

MR 

MTHTA 

MX 

MZ 

NPRT 

NULL 

OMC 

PLOSS 

PRS(I,K) 

PSI 

R(I,K) 


l/(T-l), eq. (5) 
temporary storage 
T (input variable) 

increment along quasi -orthogonal in r-direction 

increment along quasi -orthogonal in z-direction 

subscript to indicate number of quasi-orthogonal, 1 at inlet and 
MX at outlet 

code number for use by subroutine CORTOT 

number of iterations to be performed after ERROR is less than TOLER 
or after ERROR has started to increase (input variable); if 
ITER = 0, data will be printed for every iteration; if ITER > 0, 
data will be printed only for final iteration 

subscript 

subscript to Indicate number of streamline, 1 at hub and KMX at 
shroud 

number of streamlines (input variable) 

KMX - 1 

number of r values of TN in thickness table (input variable) 
number of values of THTA in table of 9 against z (input variable) 
number of fixed lines (input variable) 

number of z values of OT in thickness table (input variable) 
data is listed for every (NPRT^h streamline (input variable) 
dummy variable, not used 
1. - CORFAC 

Ap" at outlet (input variable), eq. (5) 

P 

V, eq. (7) 
r 
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(input variable) , eq. (l4) 


RB 

RC 

RH(I) 

RHO 

ROOT 

RS(I) 

RUWO 

SA(I,K) 

SAL(I,K) 

SB(I,K) 

SBETA(l,K) 

SC(I,K) 

SD(I,K) 

SFACT 

SM( I,K) 

SRW 

T 

TEMP 

THTA(J) 

TN(J,K) 

TOLER 

TP 

TPP1P 


r b 

Vr c 

r -coordinate of hub (input variable) 
p_!g (input variable) 

V2 

r -coordinate of shroud (input variable) 

integer , run number 

k, eq. (2) 

sin a 

C, eq. (2) 

sin 0 

B, eq. (2) 

D* sq* (2) 

blade multiplier to allow for splitter blades (input variable) 
distance from inlet along meridional streamline 

integer (input variable) that will cause subroutines to write out 
data for certain values, used in debugging; SRW =* 13 causes 
SPLINE to write 

t n ( interpolated) 

(input variable) 

0 (as function of z) (input variable), blade shape (fig. 6) 

t n (input variable), first subscript refers to z-coordinate, 
second subscript refers to r-coordinate 

if maximum calculated streamline correction is less than TOLER, 
iterations are considered to have converged and desired output is 
printed (input variable) 

r S0/S z 

T"/TP eq. (4) 
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TQ 

tt( i,k) 

TYPE 


TIP 

W 

WA( I,K) 
WAS 
WASS 
WT 

WTFL(K) 

WTHRU 

WTOLER 

WTR(I,K) 

XN 

xr(j) 

xt(j) 

xz(j) 

Z(I,K) 

ZH(I) 

ZS(l) 

Z SPLIT 


r S0/Sr 
t Q , eq. (9) 

integer (input variable), used as code to indicate how arrays DR, 

WA, Z, and R are given initially 

0 - These quantities will be calculated by program 

1 - These quantities are given as input on binary cards 

2 - Quantities just computed for previous case will be used for 

next case (Used only when more than one case is calcu- 
lated on single computer run) 

T/T l, eq. (3) 

o) (input variable) 

W, eqs. (l) and (13) 

W*, eq. (13) 

W**, eq. (13) 

total weight flow 

calculated total weight flow between hub and K"^ streamline, 
eq. (6) 

W n , eq. (7) 

If | WTFL(KMX) - WT| < WTOLER (input variable) , then velocity dis- 
tributions used for computing eq. (6) is accepted as solution to 
eq. (l) 

W t , eq. (10) 

N (input variable) 

r -coordinate of TN in thickness table (input variable) 
z-coordinate of THTA for blade shape (input variable) 
z-coordinate of TN in thickness table (input variable) 
z 

z-coordinate of hub (input variable) 
z-coordinate of shroud (input variable) 
z-coordinate where splitter ends (input variable) 
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Fortran Program Listing 


COMMON SOW 

DIMENSION AL(21,21 ) ,8FTA(21 .21 ) ,CAL(21 .21 ) ,rBFTA(21 .21 ) . 

1 CUR V ( 21. 21). DN (21*21). PPS (21.21). R<21. 21). Z(21. 21). SM (21.21). 

2SA( 21 ,21 ) ,SB( 21 ,21 ) ,SC( 21 ,21 ) ,SD( 21 ,21 ) *SAL( 21.21 ) .SBFTA ( 21 .21 ) . 
3TN< 21.21 ) .TT ( 21 ,21 ) ,WA( 21 ,21 ) ,WTR ( 21 .21 ) 

DIMENSION AB(21) »AC(21).AD(21) »3A(21) .DELBTA ( 21 ) .DRDM ( 21 ) . 

1DTDR (21).DTDZ(21) . DWMDM ( 2 1 ) ,DWTDM(21) . RH ( 2 1 ) . R S ( 2 1 ) , ZH ( 2 1 ) ,ZS(21) . 
2THTA( 21 ) ,WTFL<21 ) »XR (21 ) ,XT (21 ) ,XZ (21 ) 

INTEGER RUNO.TYRF.BCDP.SPW 
RUN 0=0 

10 READ ( 5, 1010)MX,KMX,MP,MZ*W,WT,XN,GA M »AP 
ITNO = 1 
RUN0=RUN0+1 
WRITE (6,1020) RUNO 

WRITE ( 6 ♦ 1010 )MX .KMX .MR ,MZ .W. WT ,XN .GAM , AR 

READ (5,1010) TYPE ,BCDP , SR W, NULL , TEMP, ALM,RHO , TOLER ,PLOSS ,WTOLER 
WRI TE ( 6 , 1010 ) TYPE ,BCDP , SR W, NULL , TEMP ,ALM,RHO , TOLER .PLOSS » WTOLER 
READ ( 5,1010)MTHTA,NPRT,ITER,NULL,SFACT,ZSPLIT,BETIN,PB,C0RFAC 
WR I TE ( 6 , 1 010 ) MTHT A , NPRT ,ITER,NULL,SFACT,ZSPLIT,BETIN,RB,COREAC 
READ( 5,1030) (ZS( I ) ,1=1 ,MX) 

WRITE(6.1030) ( Z S ( I ) ,1 = 1 »MX) 

READ( 5,1030) (ZH(I ) ♦ I = 1 .MX ) 

WRITE( 6,1030) (ZH( I ) .1=1 ,MX) 

READ ( 5,1030) (RS( I ) ,1=1 ,MX) 

WRITE(6,1030)(PS(I),I=1,MX) 

READ! 5 , 103 0) ( RH ( I ) , I = 1 , M X ) 

WRITE (6, 1030) (RH( I ) ,1=1 .MX) 

DO 20 1=1, MX 
ZS( I )=ZS( I ) / 1 2 , 

Z H ( I ) = ZH ( I ) / 1 2 . 

RS( I ) = RS ( I ) / 1 2 • 

20 RH< I ) = PH ( I ) / 1 2 , 

I F ( TYPE.NE.O) GO TO AO 

WA (1,1) = WT/RHO/ ( ZS( 1 )-ZH( 1 ) ) /3. 1A/(RS(1 )+RH( 1 ) ) 

DO 30 1=1, MX 

RN( I »KMX)=SQRT ( ( ZS ( I ) -ZH ( I ) )**2+(RS( I )-RH( I ) ) **2 ) 

DO 30 K = 1 , KMX 

DN( I ,K ) =ELOAT (K-l ) /FLOAT ( KMX-1 ) *DN ( I ,KMX ) 

WA ( I , K ) =WA ( 1,1) 

Z(I,K)=DN(I,K)/DN(I,KMX)*(ZS(I)-ZH(I) )+ZH( I ) 

30 R(I»K)=DN(I,K)/DN(I,KMX)*(RS( I) -RH ( I ) )+RH( I ) 

GO TO 50 

AC IF( TYPF.NE. 1 ) GO TO 1A5 

CALL BCREADfDN(l.l) »DN ( 2 1 ,21 ) ) 

CALL BCREAD ( WA ( 1 , 1 ) , WA ( 2 1 , 2 1 ) ) 

CALL BCREAD ( Z { 1 , 1 ) , Z ( 2 1 , 2 1 ) ) 

CALL BCREAD ( R ( 1 » 1 ) » R ( 2 1 » 2 1 ) ) 

WRITE ( 6 , 1 OAO ) 

50 READ ( 5,1030) (THTA( I ), 1=1 .MTHTA) 

WRITE (6,1030) (THTA( I ) , 1 = 1, MTHT A) 

READ ( 5 , 1030 ) ( XT ( I ) , I =1 , MTHTA ) 
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non n n n 


WRITE! 6, 1030) (XT (I ) , I = 1 >MTHT A ) 

DO 60 K=1,MR 

READ (5,1030) (TN(I ,K) ,1 = 1, MZ) 

60 WRITE (6,1030) (TNI I ,K) , 1 = 1, MZ ) 

READ (5,1030) (XZ(I), 1 = 1, MZ) 

WRITE (6,1030) (XZ(I ) ,1=1, MZ ) 

READ ( 5,1030) (XR( I ) , 1=1 ,MR ) 

WRITE (6,1030) (XR(I ) ,1 = 1, MR) 

C 

C END OF INPUT STATEMENTS 
C 

C SCALING-CHANGE INCHES TO FEET AND PS I TO LB/SQ FT, 
C CALCULATE CONSTANTS 

C 

70 DD 90 K=1,MR 
DO 80 1=1, MZ 

80 TN ( I , K ) = TN( I , K ) / 1 2 • 

90 XR(K) = XR ( K ) / 1 2 • 

DO 100 1 = 1, MZ 
IOC XZ( I ) = XZ ( I ) / 1 2 • 

DO 110 K= 1 , KMX 
110 SM ( 1 , K ) = 0 . 

BA ( 1 )=0. 

DO 120 K = 2 , K M X 

12C BA(K) = FLOAT ( K - 1 )*WT/FLOAT (KMX -1) 

DO 130 1 = 1, MX 

130 DN( I , 1 ) =0. 

DD 1 AO I = 1 >MTHTA 
140 XT ( I )=XT( I ) / 1 2 . 

ROOT = SORT (2,0) 

145 CONTINUE 

TOLER =T0LER/12. 

RB=RB/12. 

ZSPLIT = Z SPLIT/12. 

PLOSS=PLOSS*144. 

Cl = SORT (GAM*AR*TEMP) 

WRITE (6,1050) Cl 
KMXM1 = KMX-1 
CP= \R*GAM/ (GAM-1 . ) 

EXPON = 1. /(GAM-1.) 

BETIN = -BETIN/57. 29577 
RINLET = (RSC 1 )+RH(l ) )/2. 

CEF=S IN (BETIN) /COS (BETIN) /R I NLET / < R I NLET-RB ) * *2 
ERR OR = 100 000 • 

BEGINNING OF LOP D FOR ITERATIONS 

150 I F ( ITER.EQ.O) WRITE (6,1060) ITNO 
I F ( ITER.EQ.O) WRITE (6,1070) 

ERR0R1=ERR0R 
ERR0R=0. 

START CALCULATION OF PARAMETERS 
DO 230 K = 1 , K M x 


INITIALIZE, 
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n n n n 


DO 160 1 = 1, MX 

AB ( I ) = < Z < I ,K )-R ( I ,K ) ) /ROOT 
160 AC ( I) = ( Z ( I ,K )+R ( I ,K ) ) /ROOT 

CALL S D L INF ( AB , AC , MX , AL ( 1 , K ) , CUR V ( 1 , K ) ) 

DO 170 1 = 1, mx 

CUPV (1,0 = CURV ( I ,K ) / { 1 ,+AL ( I , O **2 )**1 .5 
AL (I*<> = ATAN(AL( I ,0 )-. 785398 

CAL (I,<) = COS ( AL ( I , K ) ) 

170 SAL (1,0 = S I N ( AL (1,0) 

DO 180 1=2, MX 

180 SM ( 1,0 = 5M( I-1»K)+SQRT ( (Z ( I ,K)-Z(I-1»0 ) **2 + ( R ( I » O - P ( 1-1,0 ) ** 

1 2 ) 

190 CALL SPLDFR ( XT ( 1 ) ,THTA ( 1 ) ,MTHTA,Z ( 1 ,K ) ,MX ,OTD7 ( 1 ) ) 

DO 220 1=1, mx 

CALL L I N I NT ( Z ( I , K ) ,R( 1,0 ,XZ,XR,TN, 21,21, T) 

I F ( P ( I ,<) . LE « RB ) GO TO 2 00 
DTDR ( I )=CEF*(R(I,0~RB)**2 
GO TO 210 
200 DTDR ( I ) =0. 

210 TQ = R (1,0 *DTDR ( I ) 

TP = R ( I , O *DTDZ ( I ) 

TT ( I » O = T*SQRT ( l.+TP*TP) 

BETA(I,0=ATAN(TP*CAL(I »O+T0*SAL { I ,K ) J 
SBETA ( I , O = SI N ( BFTA ( 1,0) 

CBET A ( I , K ) = COS ( BETA (1,0) 

SA( I ,K)=CBETA( I ,K) **2*CAL (1,0 *CURV( 1,0 -SBETA (1,0 **2/R ( I ,K ) + 

1 S AL ( I »0*CBETA( I »0*SBETA( I ,K ) *DTDR ( I ) 

SC ( I , K ) =-S AL ( I , K ) *CBET A ( I ,0**2*CURV( I ,0+SAL( I ,0 *CBFTA ( I ,0 
1 *SB ETA ( I ,K)*DTDZ ( I ) 

AB ( I ) =WA ( I , K ) *CB E TA ( I ,K ) 

220 AC ( I) =WA( I ,K ) *SBETA ( I ,K ) 

CALL S D L I NE ( SM ( 1 , K ) , AB , M X , DWMDM , AD ) 

CALL SPL I NE ( SM ( 1,0 , AC ,MX , DWT DM , AD ) 

I F( ( ITER.LE.O) .AND. (MOD (K-l ,NPRT ) .EQ.O) ) WRITE (6,1080) K 
DO 230 1=1, MX 

SB ( I * O = SAL ( I ,K)*CBETA( I ,K)*DWMDM( I ) -2 . *W*SBETA ( I , K ) +DT DR ( I ) * 

1R( I ,0*CBFTA( I ,< ) * ( DWTDM ( I )+2.*W*SAL( 1,0) 

SD ( 1,0 =CAL (1,0 *CBETA ( I , K) *DWMDM( I J+DTDZ ( I ) * 

1 R ( I , K ) *CBETA( I ,K ) * ( DWTDM ( I ) +2.*W*SAL< I ,K ) ) 

IF( ( ITER.GT.O) .OR. ( MOD < K- 1 , NPR T ) .NF.O) )G0 TO 230 
A= AL ( I ,K)*57. 29577 
B= SM (1,0*12, 

E= TT ( I , K ) * 1 2 , 

G= BETA ( I ,0*57.29577 

WRITE (6,1090) A,CURV( I ,K) ,B,G,E, SA ( I , O , SB ( I ♦ O , SC ( I , O , SD ( I , K ) 
230 CONTINUE 

END OF LOOP - DARAMETER CALCULATION 

CALCULATE BLADE SURFACE VELOCITIES (AFTER CONVERGENCE) 

I F ( ITER.NE.O) GO TO 260 
DO 250 K = 1 , KMX 

CALL SPLINE ( SM ( 1 ,K ) , TT ( 1 ,0 , MX , DELBT A , AC ) 

A = XN 
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n n n o 


DO 240 1 = 1, MX 

24 J AB( I ) = { R ( I ,0*W+WA( I »0*SBETA( I, 0)*(6. 283186 *P <1,0/ A-TT( I.K) ) 
CAL'. SPLINE ( SM ( 1 , O , AB ,MX ,DRDM, AC ) 

IF (SFACT.LF. 1.0) GO TO 245 
A = SF ACT *XN 
DO 244 1 = 1 , MX 

2 44 AB ( I ) = ( R ( 1,0 * W + W A (1,0 *SRFT A ( I , K ) ) * ( 6 . 283 1 86 *P ( I , K ) / A-TT ( I ,0 ) 
CALL SPLINE ( SM ( 1 ,0 , AB , M X , AD ,AC) 

245 DO 250 1=1, MX 

BETAD = BETA ( I ,K) -DFLPT A ( T ) / 2 • 

BETAT = BET AD + DELBT All) 

COSBD = CO S ( BET AD ) 

COSBT = COS ( BET AT ) 

IF(Z( I , O .LT.Z SPLIT) DRDM(I) = AD { I ) 

WTP ( I »O=C0S3D*CQSBT/( COSBD+COSBT ) *< 2.*WA ( I ,K ) /COSBD+P (1,0 *W* 

1 (BETAD-BFTAT)/CRETA( I »O**2+0RD M ( I ) ) 

250 CONTINUE 

END OF BLADE SURFACE VELOCITY CALCULATIONS 
START CALCULATION OF WEIGHT FLOW VS. DISTANCE FROM HUB 

260 DO 070 1 = 1, MX 
I ND= 1 

DO 270 <= 1 , KMX 
270 AC ( K ) = DN (1,0 
GO TO 290 

280 W A ( 1,1 ) =.5*WA( 1,1) 

293 DO 300 K = 2 * KMX 
J = K-1 

HR = R ( I , O -R ( I , J ) 

HZ = Z ( I ,K) -7 ( I , J) 

WAS = WA ( I , J )* ( 1 .+SA( I ,J )*HR + SC ( I * J ) *HZ ) +SB ( I , J )*HP+SD( I * J ) *HZ 
WASS = WA ( I , J ) +WAS* ( SA( I , K ) *HR + SC( I , K ) *HZ ) +SB ( I ,0 *HR+SD( I ,0 *HZ 
300 WA ( I , K ) = ( WAS+WASS ) / 2 . 

310 DO 340 K= 1 » KMX 

T 1 P = 1 . - ( W A ( I ,K) **2+2.*W*ALM- ( W*R (1,0 )**2 ) /2 ./CP /TEMP 
I F ( T1P.LT..0) GO TO 280 

TPP1P= 1. ( 2 .*W*ALM- ( W*R ( I ,0 )**2 ) /2 ./CP/TEMP 

DENSTY = TlP**EXPON*RHO-(TlP/TPPl D )*«-EXPnN#PLOSS/AP/TPPlP/TEMD 
1 *32. 17*SM( I ,K)/SM(MX,K) 

PRS ( I , O = DENST Y*AR*T 1 D *TEMP /3 2 . 17/144. 

I F ( Z S ( I ).LE.ZH(I ) ) GO TO 320 

PS I = AT AN ( (RS( I )-RH< I ) ) / (ZS( I ) -ZH ( I )) 1-1.5708 
GO TO 330 

320 PSI=ATAN( (ZH( I )-ZS( I ) )/(RS( I ) -RH ( I ) ) ) 

3 30 WTHRU = WA ( I ,O*C0ETA( I ,K ) *CQS ( PSI-AL ( I ,K ) ) 

A = XN 

I F ( Z ( I ,K) .LT.ZSPLIT) A = SFACT*XN 
C = 6.283186*R( I »0-A*TT ( I » O 
340 ADtO=DENSTY*WTHRU*C 

CALL INTGRKACIl ) , AD ( 1 ) ,KMX,WTFL(1) ) 

IF ( ABS( WT-WTFL ( KMX ) ) .LE.WTOLER ) GO TO 350 
CALL CONTIN ( WA ( I , 1 ) , WTFL ( KMX ) , I ND , I , WT ) 

IF (IND.NE.6) GO TO 290 
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n n n n 


350 CALL SPLINT ( WTFL , AC * KMX , BA , KMX ♦ A3 ) 

DO 360 K=1,KMX 
DELTA=ABS(AB(K)-DN( I, KM 

DN ( I » K ) = ( 1 .-CORF AC ) *DN ( I »K)+CORFAC*AB(K) 

360 IF ( DEL T A. GT. ERROR ) ERROR = DEL T A 
370 CONTINUE 

END OF LOOP - WEIGHT FLOW CALCULATION 

CALCULATE STREAMLINE COORDIMATFS FOR NFXT I TFRAT I ON 


DO OPQ K=2,KMXM1 
DO 380 1=1, MX 

Z( I ,K)=DN( I »K)/DN( I »KMX ) * ( ZS ( I )-ZH ( I ) )+ZH( I ) 

380 R(I»K)=DN(I,K) /DN ( I »KMX ) * ( RS ( I ) -RH ( I ) )+RH( I ) 

I F ( (ERROR.GE. ERROR 1 ) .OR. ( ER ROR . LE . TOLER ) ) ITER= ITER-1 
I F { ITFR.GT.O) GO TO 410 
WRITE (6,1100) 

DO 400 K = 1 »KMX,NPRT 
WRITE (6,1080) K 
DO 390 1=1, MX 

AB ( I ) = ( Z ( I , K ) -R ( I , K ) J/ROOT 
390 AC( I ) = (Z( I ,K)+P( I ,K) J/ROOT 

CALL SPLINE ( AB , AC , MX , AL ( 1 , K ) ,CUR V ( 1 , K ) ) 

DO 400 1 = 1 , MX 

CUR V ( I , K ) =CURV (I,K)/(1.+AL(I,K)**2)**1.5 
A=DN ( I ,K) *12. 

B= Z ( I ,K ) *12. 

D= R ( I ,K ) *12. 

400 WRITE (6,1110) A , B , D , W A ( I , K ) , PRS ( I , K ) , WTR ( I , K ) , CUR V ( I , K ) 
WRITE (6,1130) 

410 A = ERRRR * 12. 

WRITE (6,1120) I T NO , A 
I TNO=I TNO+1 

IF ( ITFR.GE.O) GO TO 150 

IF(BCDP.NE.l) GO TO 10 

CALL BCDUMP ( DN ( 1 , 1 ) , DN ( 2 1 , 2 1 ) ) 

CALL BCDUMP ( WA ( 1 , 1 ) , W A ( 2 1 , 2 1 ) ) 

CALL BCDUMP ( Z(l,l), Z(21,21)) 

CALL BCDUMP ( R(l,l)» R(21,21)) 


42° GO TO 10 
1010 FORMAT ( 4 I 5 ♦ 6F1 0 • 4 ) 


1020 FORMAT 
1 03 C FORMAT 
1040 FORMAT 
1050 FORMAT 
1060 FORMAT 
1070 FORMAT 
1X5HSB 
1080 FORMAT 
1090 FORMAT 
1100 FORMAT 


(8H1RUN NO.I3»10X,25HINPUT DATA CARD LISTING ) 

(7F10.4) 

( 1 0X24HBCD CARDS FOR DN,WA,Z,R ) 

( 36HK STAG. SPEED OF SOUND AT INLET = ,F9.2) 

( ///5X13HI TERATION NO. 13) 

(1H 6X5HAL 9X5HRC 9X5HSM 9X5HBET A RX^HTT 9X5HSA 
9X5HSC 9XFHSD ) 


(2X10HSTREAMLINFI3) 

(9F14.6) 

(1HL9X5HDN 15X5HZ 15X5HR 
1TR14X3HRC ) 

1110 FORMAT (6F19.6,F18.6 ) 

1120 FORMAT (18H ITERATION NO. 

1 F 1 0 . 6 ) 

1130 POP MAT (1HJ) 

FND 


15X5HWA 


15X5HPRS 14X3HW 


I 3 , 1 OX , 24HMAX • STREAMLINE CHANGE = 
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Description of Subroutines 

The subroutines SPLINE, SPLINT, SPLDER, and INTGRL are based on the spline 
fit curve (see appendix C) . SPLINE gives the first and second derivatives, 
SPLINT is used for interpolation, SPLDER is used for interpolated values of the 
derivative, and INTGRL is used for numerical integration of a function given at 
unequally spaced points. The calling sequences for these subroutines are as 
follows : 


CALL SPLINE (X,Y,N, SLOPE, EM) 


where 

X input array 

Y input array, function of X 

N input, number of X and Y values given 

SLOPE output array, first derivative, dY/dX 

EM output array, second derivative, d^Y/dX^ 


CALL SPLINT (X,Y,N,Z,MAX,YINT) 


where 

X input array 

Y input array, function of X 

N input, number of X and Y values given 

Z input array, values at which interpolated function values are desired 

MAX input, number of Z values given 

YINT output array, interpolated values 


CALL SPLDER (X, Y,N,Z,MAX,DYEK) 


where 



X 


input array 


Y input array, function of X 

N input, number of X and Y values given 

Z input array, values at which the derivative is desired 

MAX input, number of Z values given 

DYEX output array, derivatives at each Z 


CALL INTGRL (X,Y,N,SUM) 

where 

X input array 

Y input array, function of X 

N input, number of X and Y values given 

/•X(I) 

SUM(l) output array, / Y DX 

•'x(i) 

The subroutines SPLINE, SPLINT, SFLDER, and INTGRL are as follows: 
SUBROUTINE SPLINE ( X *Y *N ♦SLOPE , EM )' 

DIMENSION X( 50) ♦ Y ( 5 0 ) ♦$(50) ♦ A ( 5 0 ) ♦ B ( 5 n ) ♦ C ( 5 P ) ♦ E ( 5 0 ) ♦ W ( 5 0 ) * S B ( E> 0 ) ♦ 
1 G ( 5 0 ) , EM (50) * SLOPE ( 50 ) 

COMMON Q 
INTEGER 0 
DO 10 I = 2 ♦ N 
10 S( I )=X ( I ) —X ( 1-1 ) 

N0=N-1 

DO 20 1=2 ♦NO 
A(I)=S(I)/6. 

B( I )={S(I)+S( I+l ) )/3. 

C ( I ) = S ( I + 1 ) / 6 • 

20 Ft I ) = ( Y( I + l )-Y( I ) )/S< I + l )-(Y( I ) -Y < 1-1 ) )/S< I ) 

A ( N ) =-.5 
B ( 1 ) =1 . 

B ( N ) = 1 . 

C( 1 ) =- . 5 
F ( 1 ) =0. 

F ( N ) =0 . 

W( 1 )=B ( 1 ) 

SB ( 1 ) =C ( 1 ) /W ( 1 ) 

G ( 1 ) = C . 

DO 30 I =2 ♦N 

W( I ) =B( I )-A( I )*SB ( 1-1 ) 

SBC I )=CC I )/WC I > 


33 



30 G( I ) = ( F ( I ) — A ( I )# 0(1-1 ) ) / W ( I ) 

EM ( N ) = G ( N ) 

DO 40 1=2, N 
K = N+ 1 - I 

40 EM( K )=G(K)-SB(K)*EM(K+1 ) 

SL0PE(1)=-S(2)/6.*(2.*FM(1)+EM(2) ) + ( Y ( 2 ) -Y { 1 ) )/S(2) 

DO 5 0 1=2, N 

50 SLOPE <I)=S(I)/6.*( 2 .*EM ( I )+EM (I-1)) + (Y(I)-Y(I-1))/S(I) 

IF ( Q. EQ» 1 3 ) WRITE (6,100) N*(X(I)»Y(I)»SL0PE(I)»EM(I)»I=1»21) 

100 FORMAT (2X15HN0. OF POINTS =I3/10X5HX 15X5HY 1 5 X5 HSL0PE1 5X 5H 
1EM /(4F20.P)) 

RETURN 

END 


SUBROUTINE S°LINT ( X , Y , N , 7 , M AX , Y I NT ) 

DIMENSION X(50),Y(5P),S(50),A(50),B(50),C(50),F(50),W(50),SB(50), 
1G ( 5 0 ) , EM I 50) ,Z ( 50 ) ,YINT(5 0) 

COMMON Q 
INTEGER 0 
DO 10 1=2, N 
10 S< I) =X ( I ) -X ( I -1 ) 

N 0 = N - 1 

DO 20 1=2, NO 
A t I ) =S ( I ) / 6 • 0 
B ( I ) = (S( I )+S( I + 1 ) 1/3.0 
C ( I ) =S ( I +1 ) / 6 . 0 

20 F( I ) = (Y( I +1 ) — Y ( I ) )/S( I + 1 )-(Y( I ) — Y ( I— 1) )/S( I ) 

A (N ) =-. 5 
B ( 1 ) =1 .0 
B ( N ) = 1 . 0 
C(1 ) =- . 5 
F ( 1 ) =0.0 
F ( N ) =0.0 
W( 1 )=B( 1 ) 

SB ( 1 ) =C ( 1 ) /W < 1 ) 

G( 1 ) =0.0 

DO 30 1=2, N 

W( I )=B( I ) —A ( I ) * S B ( 1-1 ) 

S B ( I ) = C (I ) / W ( I ) 

3 0 G( I ) = ( F( I ) —A C I )*G( 1-1 ) )/W( I ) 

EM ( N ) = G ( N ) 

DO 40 1=2, N 
K = N + 1- I 

40 EM( K ) =G( K ) -SB < K) *EM(K+1 ) 

DO 90 1=1, MAX 
K = 2 

I F ( Z ( I ) —X ( 1) ) 60,50,70 
50 Y INT ( I ) =Y ( 1 ) 

GO TO 90 

60 I F ( Z ( I ).LT.(1.1*X(1)-.1*X(2) ) ) WRITE (6*1 000 ) Z ( I ) 

GO TO 85 

1000 FORMAT (17H OUT OF RANGE Z =F10.6) 

65 I F ( Z( I ) .GT. ( 1.1*X (N)-.1*X(N-1 ) ) ) WRITE (6*1 000 ) Z ( I ) 

K=N 

GO TO 05 
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70 I F ( 7. ( I ) -X ( K ) ) 85*75*80 
75 YINTf I >=YfK) 

GO TO 90 
83 K=K+1 

IFCC-N) 70,70,65 

85 YINT(I) = EM( K-l ) * ( X ( K )-Z ( I ) ) **3/6. /S <K )+EM( K )* (Z ( I ) -X ( K-l ) )**3/6 . 
1/S< K ) + ( Y( K ) /S (K )-EM(K)*S (K ) /6 . )* ( Z ( I )-X (K-l ) ) +( Y( K-l ) /5 ( K )-FM(K-l ) 
2*S ( K ) /6. ) * ( X ( K )-Z < I )) 

90 CONTIMUF 

I F ( Q.EQ. 16 ) WR I TE (6 , 1010 ) N * M AX , < X (I ) » Y ( I ) » Z ( I ) * Y I NT ( I ) , I = 1 , N ) 

1010 FORMAT (2X21HN0. OF POINTS GIVEN =»I3»30H, NO. OF INTERPOLATED POI 
1NTS = , I 3 * / 1 OX 5HX 15X5HY 1 2X 1 1HX- I NTERPO L . 9X 1 1 H Y- I N TER POL . / ( 4 

2E20.8 ) ) 

100 RETURN 
END 


SUBROUTINE SPLDFR < X , Y ,N ,Z ,MAX ,DYDX ) 

DIMENSION X ( 5 0 ) , Y ( 5 0 ) * S ( 5 0) , A ( 50) ,8(50) ,C(50) ,F(50) ,W( 50 ) ,SB(50) , 
1G( 50) ,EM( 50) ,Z ( 50 ) ,DYDX( 50) 

DO 10 1=2, N 
10 S ( I ) =X ( I ) -X ( 1-1 ) 

N 0 = N - 1 

DO 20 1=2, NO 
A ( I )=$( I ) / 6 • 0 
B( I ) = (S( I )+S( 1+1 ) ) / 3 • 0 
C ( I ) =S< 1+1 ) /6.0 

20 F( I ) = (Y( I + 1 ) - Y { I ) )/S( 1 + 1 ) - ( Y ( I ) - Y '( 1-1 ) )/S( I ) 

A ( N ) =-. 5 
B ( 1 1 = 1*0 
B ( N ) = 1 .0 
C ( 1 ) =-. 5 
F ( 1 )=0.C 
F ( N ) = 0 • 0 
W ( 1 ) =B ( 1 ) 

SB(1)=C(1)/W(1) 

G ( 1 ) = 0 . 0 
DO 30 1=2, N 

W( I )=B( I ) -A ( I ) *SB ( 1-1) 

SB ( I 1 = C ( I ) /W ( I > 

30 G(I) = (F(I)-A(I) *G ( I - 1 ) ) / W ( I ) 

EM ( N ) = G ( N ) 

DO 40 1=2, N 
K = N + 1 - I 

40 EM(K)=G(K)-SB(X)*EM(K+1) 

DO 90 1=1, MAX 
K = 2 

IF(T(I)-X(1)) 60,70,70 

60 WRITE ( 6 , 1 000 ) 7 ( I ) 
l^OO FORMAT (17H OUT OF BLADE Z =F10.6) 

GO TO 85 

6 5 WRITE ( 6,1000 )Z( I ) 

K=N 

GO Tf 45 

70 I F ( Z ( I ) — X ( K > ) 85,85,80 
o o K = K + 1 
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n n 


IF(Y-N) 70,70,65 

8 5 DYDX ( I)=-EV(K-1)*(X(K)-Zm ) **2/2 .0/5 ( K ) + EM ( K ) * ( X ( K-l )-Z< I) >**2/2. 

10/S(K) + (Y(K) -Y ( K - 1 ) ) /S (<■)-( EM (K)-EM(K-l ) )*S(K )/&. 

90 CONTINUE 
ICO RETURN 
END 


SUBROUTINE INTC-PL (X,Y,N,SUM) 

DIMENSION X ( 5 C ) ,Y (50) »S ( 50 ) , A ( 5 0 ) * B < 5 0) ,C(50) ,F ( 50 ) »W( 50 ) ,SB ( 50 ) * 

1 G { 50 ) , EM ( 50 ) .SUM ( 50 ) 

DIMENSION X ( 50) , Y ( 5 0 ) »S ( 50) »A ( 50) ,3(50) ,C ( 50) ,F ( 50) »W( 50) , SB ( 5 0 ) , 

1 G ( 50) ,EM( 50) , SUM ( 50 ) 

DO 10 I =2 »N 
10 S( I )=X( I )-X( I — 1 ) 

N 0 = N - 1 

DO 20 1=2, NO 
A ( I ) =S ( I ) /6. 0 
B ( I ) = ( S ( I >+S( 1 + 1 ) ) / 3 . 0 
C ( I ) =S( 1+1 ) /6.0 

20 F( I ) = (Y( I +1 ) — Y ( I ) )/S( 1 + 1 ) -( Y ( I ) — Y ( I — 1 ) )/S( I ) 

A ( N ) =- . 5 
B< 1 )=1 .0 
B ( N ) =1 .0 
C ( 1 ) =-. 5 
F ( 1 ) =0.0 
F ( N ) =0.0 
W ( 1 ) = B ( 1 ) 

SB ( 1 ) =C ( 1 ) /W( 1 ) 

G( 1 ) =0.0 

DO 30 1=2, N 

W( I ) =B ( I) -A ( I ) *SB ( 1-1 ) 

SB( I )=C( I )/W( I ) 

30 G( I ) = (F( I ) — A ( I ) *G ( I — 1 ) )/W( I ) 

EM ( N ) =G ( N ) 

DO 40 I = 2 , N 
K=N+1- I 

40 EM ( r ) =r, ( K) - SB < K ) *EM { K+l ) 

SUM ( 1 ) =0.0 

DO 50 K = 2 » N 

50 SUM ( K ) = SUM ( K-l )+S ( K ) * ( Y (K )+Y( K— 1 ) ) /2 .0-S ( K ) **3* ( EM ( K ) +FM ( K-l ) ) /2 

14.0 
RETURN 
END 


The subroutine LININT performs linear interpolation of a function of two 
variables. It is used here to obtain interpolated values of normal blade 
thickness t n from a table of thickness values given as input. The calling 
sequence for LININT is as follows: 

CALL LININT (X1,Y1,X,Y,TN,MX,MY,F) 


where 
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XI input, x-coordinate of point for which interpolated function value is 
desired 

Y1 input, y-coordinate of point for which interpolated fraction value is 
desired 

X input array, x-coordinates at which function values are specified 

Y input array, y-coordinates at which function values are specified 

TN input two-dimensional array, function of x and y, first subscript 
refers to x-coordinate 

MX input, number of x values given 

MY input, number of y values given 

F output, interpolated value 

The subroutine LININT is as follows: 


SUE ROUTINE LININT(X1»Y1.X,Y,TN*MX,MV»F) 
common < 

DIMENSION X ( MX ) , Y (MY ) , TN ( MX »MY ) 

DO 10 J3= 1 »MX 

10 I F { X1.LE.X( J3 ) )G0 TO 20 
J3 = MX 

20 DO 30 J4 = 1 » w Y 
30 I F ( Y1.LE.Y(J4) )G0 TO 40 
J4 = MY 

40 J 1 = J3- 1 
J2= J4-1 

E D S 1 = ( Xl-X ( J1 ) )/(X(J3)-X(Jl) ) 

EPS2=(Y1-Y(J2) ) / ( Y ( J4 ) — Y ( J2 ) ) 

EPS 3 = 1 . —EPS 1 
F D S 4= 1 . — F PS2 

F=TN(J1,J2)*EPS3*EPS4+TN(J3*J2) *EPS 1 *EPS4+TN (J1»J4)*EPS2*EPS3+ 
1TN( J3,J4)#EPS1*EPS2 

I F ( K • E0 • 1 4 ) WRITE(6*1)X1*Y1»F*J1*J2»EPS1 ,EPS2 
1 FORMAT ( 8 H L I N I NT 3 F 1 0 . 5 , 2 I 3 * 2 F 1 0 . 5 ) 

F=0 

RETURN 

END 


The subroutine C0MTIN is used to predict the hub velocity to be used in 
the next iteration to satisfy continuity of flow (eq. (6)) between hub and 
shroud. An initial estimate is furnished by the main program, say Wp (see 
fig* 12) . C0MTIN furnishes the next estimate W2 by linear interpolation or 
extrapolation from the origin. Subsequent estimates are obtained by linear 
interpolation from the two previous estimates. 
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Figure 12. - Method used by subroutine CONTIN to determine relative hub velocity. 


If there is choked flow, 
there is no solution of equa- 
tion (l) that will also satisfy 
continuity (eq. (6)). In this 
case, CONTIN will find the hub 
velocity that gives the maximum 
calculated weight flow. It 
should be noted that CONTIN 
does not calculate the weight 
flow; this Is calculated by the 
main program. CONTIN stores 
information from up to three 
previous iterations to assist 
in predicting the next value 
to be used for the hub velocity. 
The calling sequence for CONTIN 
Is as follows : 


CALL CONTIN (WA,WTFL,IND,I,WT) 


where 

WA input and output; as input, hub relative velocity used to calculate 
latest weight flow and as output, velocity used for next iteration 

WTFL input, calculated weight flow based on input value of WA 

IND input and output; main program sets IND * 1 to indicate start of weight- 
flow calculation for new quasi -orthogonal and CONTIN changes value of 
INI) for following iterations to indicate procedure followed in calcu- 
lating new hub velocity 

I input, number of quasi-orthogonal used by subroutine CONTIN in WRITE 

statement if there is choked flow 

WT input, total weight flow 

The subroutine CONTIN is as follows: 

SUBROUTINE CONTIN ( WA » WTFL » I NO » I » WT ) 

DIMENSION SD p ED ( 3 ) .WEIGHT (3 1 
135 GO TO ( 1 4 0 * 1 5 0 * ? 1 C » ? 7 o , 3 ?0 ) » I ND 
l&c s D FFP ( 1 ) = WA 

WEIGHT! 1 ) = WTFL 
WA = WT/WTFL*WA 
IND = 2 
RETURN 

150 IF ( (WTFL- WEIGHT ( 1) )/(WA-SPEED(l ) ) ) 180*180»16C 

160 SPFFDI 2 ) = WA 

WA = ( WT-WTFL ) / ( WTFL- WE I GHT (11) 

1 * ( WA-SPFFD ( 1 ) ) +WA 
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IF ( ARS(WA-SPFFp(2) )-l A 0.01 166,166,161 

161 IF(WA-SPEED(2> >163,163,162 

162 WA = S°EED( 2 ) +100.0 
GO TO 166 

16 3 WA = SPEED ( 2 ) -100.0 
166 SPEED(l) = S D EFD ( 2 ) 

WEIGHT! 1) = WT C L 
RETURN 

170 WRITE ( 6 , 1000 ) I , W T F L 
IND = 6 
RETURN 
180 IND = 3 

IF ! WTFL.GE.WT ) GO TO 14 
IF ( SPFFD ( 1 ) -WA ) 190,200,200 

1 9 P S D F F D ( 2 ) = SPFFD ( 1 ) 

SPEED ( 1 ) = 2 . C*SPEED ( 1 ) -WA 

SPEED! 3) = WA 

WEI GHT ( 2 ) = WEIGHT! 1) 

WEIGHT! 3) = WTFL 
WA = SPEED ( 1 ) 

RETURN 

200 S D E E D ( 2) = VIA 

SPEED! 3) = SPEFD(l) 

SPFFD! 1) = 2.0*WA-SPFED( 1 ) 

WEIGHT (2) = WTFL 
WEIGHT (3) = WEIGHT! 1) 

WA = SPEED ( 1 ) 

RETURN 

210 WEIGHT! 1) = WTfl 

IF (WTFL.GF.WT) GO TQ 14 
IF (WEIGHT ( 1 )-WEIGHT<2 ) ) 2 30,380,220 
220 WEIGHT!?) = WEIGHT! 2) 

W E I G H T ( 2 ) = WEI GHT ( 1 ) 

SPFFD ( 3 ) = SPEFD ( 2 ) 

SPEED ( 2 ) = SPEED(l) 

SPEED(I) = 2.0* SPEED! 2 ) -SPEED (3 ) 

WA = SPEFD! 1 ) 

RETURN 

23f IF ( SPFFD ( 3 )-SOFED! 1 ) -10.0) 170,170,240 
240 IND = 4 

245 IF ( WF I GHT ( 3 ) -WEIGHT ( 1 ) ) 260,260,250 
250 WA = ( SPEED! 1 J+SPEED ( 2 ) ) /2.0 
RETURN 

260 WA = ( SPEED! 3 )+SPEED( 2 ) ) /2.0 
RETURN 

270 IF ( SPEED(3 ) -SPEED! 1 )-10.0) 170,170,280 
280 IF ( WTFL- WE I GHT ( 2 ) ) 320,350,290 
290 IF (WA-SPEED12 ) ) 310,300,300 
300 SPEED! 1 ) = SPEFD ( 2 ) 

SPEED! 2) = WA 

WE I GHT ( 1 ) = WE I GHT ( 2 ) 

WE I GHT ( 2 ) = WTFL 
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GO T n 245 

310 SPEED ( 3 ) = S D FFD ( 2 ) 

C PFFD(?) = WA 
WEIGHK3) = WFIGHT(2) 

WEI GHT ( 2 ) = WTFL 

GO TO 245 

320 IF ( WA-SPFED ( 2 ) ) 340,330,330 

330 WEIGHT! 3) = WTFL 
SPEED! 3) = WA 

GO TO 245 

Q 4 P WEIGHT! 1 ) = WTFL 
SPEED! 1) = WA 
GO TO 245 
350 IND = 5 

IF (WA-SPFED! 2 ) ) 380,360*^60 

360 SPEED! 1 ) = SPEED ( 2 ) 

WEIGHT ( 1 ) = WEIGHT ( 2 ) 

SPEED!?) = ( corrn ( ] J+SDFCP ( ? ) ) /? , 

WA = S n F F D ( ? ) 

RETURN 
t 7 C IND = 4 

WEIGHT! 2) = WTFL 
WA = ( SPEED! 1 )+S D EED( 2 ) ) /2.0 

RETURN 
3 A 3 IND = 5 

3 FO WEIGHT! 3) = WEIGHT! 2) 

S D E FD ( 3 ) = S PE RD ( 2 ) 

S D F r D ( 2 ) = ( SPFGD! 1 )+S°EFD( 3 ) ) /2. 

WA = S D EFD ( 2 ) 

RETURN 

1P00 FORMAT ( / 1 2H FIXED LINE I2,12H, MAX WT = E10.6) 
END 


Sample Output from Program 

The output given here is the listing for the case used in the numerical 
example. It will he noted that there is an exact listing of all input data 
cards at the beginning of the listing. This is followed by the maximum calcu- 
lated streamline change for each iteration, which is used as the criterion for 
convergence. After 47 iterations, there is convergence within the specified 
limit of 0.001-inch maximum streamline change. At this time, streamline co- 
ordinates are printed together with the velocity and pressure at each point. 
This is followed by another iteration to give additional information of inter- 
est, such as a, (3, and the parameters A, B, C, and D from equation (2). 
Since it indicates the smoothness of the streamline at a glance, the streamline 
curvature is also printed out. The velocities and the pressures are computed 
again on the final iteration so that the variation of these quantities on the 
final iteration can be checked. 
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RUN NJ. 1 
10 11 
0 0 

13 2 

0.3000 
0.6130 
0. 

0. 7300 
2 .2500 
1.4765 
2.2500 
0.6600 
0 . 

-0.0517 
-0.1000 
0.6000 
0.3850 
0.1750 
0.3750 
0.1650 
0.3650 
0.1550 
0.3550 
0.1460 
0.3450 
0.1360 
0.3330 
0.1270 
0.3200 
0.1170 
0.3100 
0.1070 
0.3000 
0.0980 
0.2800 
0.0880 
0.2510 
0.0790 
0.2230 
0.0700 
0.1940 
0.0600 
0.1 660 
0.0500 
0.1370 
0.0400 
0.1090 
0.0300 

o.oaoo 
0.0200 
- 0. 1000 
0.600 0 
0.6500 

1 .3500 

2 .0500 


INPUT DATA CARD LISTING 


1 7 13 

5390.0000 

0.9840 

-0 0 

592.0000 

155.3000 

2 I 

1 .0000 

- 1 .0000 

0.3400 

0.3948 

0.4810 

0.9100 

1.0000 


-0.0270 

-0.0530 

-0.0540 

0.9100 

1.0400 


2.0520 

1.8610 

1 .6800 

1*4751 

1.4750 


2.0150 

1.7630 

1.4120 

0.6750 

0.6750 


0. 

0. 

-0.0004 

-0 .0972 

-0.1632 

-0.2487 

0. 

0.1000 

0.2000 

0.7000 

0.8000 

C . 9000 

0.3580 

0.3220 

0.2900 

0.1500 

0.1280 

0.1080 

0.3450 

0.3110 

0.2800 

0. L 400 

0. 1190 

0.1000 

0.3330 

0.3000 

0.2 700 

0.1310 

0 . i 100 

0.09 20 

0 . 3 2 1 1 j 

0.2900 

0.2590 

0.1210 

0.1010 

0 . 0850 

0.3100 

0.2800 

0.2490 

0.1130 

0.0930 

0.0790 

0.3000 

0.2690 

0.2380 

0.1040 

0.0860 

0.0730 

0.2900 

0.2580 

0.2270 

0.0950 

0.0780 

0.0680 

0.2790 

0.2470 

0.2 150 

0.0870 

0.0730 

0 . 0620 

0.2680 

0.2360 

0.2040 

0.0790 

0.068 C 

0.0580 

0.2570 

0.2240 

0.1930 

0.0730 

0.0620 

0. C 540 

0.2460 

0.2130 

C . 1 820 

0.0670 

0.0570 

0.0500 

0.2230 

C • 2020 

C . 1 700 

0.0610 

0.0520 

U . U 4 oO 

0 . 1940 

0. 1900 

C . 1 590 

-0. 

-0. 

-0 • 

0 * 1 6 oO 

0. 1660 

0. 1480 

-0. 

-0. 

-0. 

0.1370 

0.1370 

0.1370 

-0 • 

-0. 

-0. 

0. 1090 

0.1090 

0. 1 050 

-G. 

-0. 

-o. 

0.0800 

0.0800 

0.0800 

-0. 

-0. 

-o. 

0. 

0.1000 

0.2000 

0.7000 

0.8000 

0. 9000 

0.7500 

0.8500 

0.9500 

1.4500 

1.5500 

l .6500 

2. 1500 

2.2500 



13.0000 

l .4000 

1715.0000 

0. 1941 

0.001 0 

2 . 5000 

-35.0000 

1 . 7500 

0. ID JO 

0.5412 

0.6 193 

0. 70 «0 

0.0 09 0 

0.1370 

0.4100 

1.6000 

l .5 34 1 

1.4960 

1 .2080 

1.0010 

U . 7 7 90 

-0.0027 

-0.0090 

-0.0240 

-0. 3512 

-0.4660 


0. 3000 

0.4009 

0.5000 

1.0000 

1 . 1000 


0.2600 

9.2300 

0. 20 10 

0.0920 

0. 0/90 


0.2500 

0.2200 

0. 19 10 

0.0840 

0.0740 


0.2400 

0.2100 

0. la 10 

0.0780 

0.0690 


0.2290 

0.2000 

0.1710 

0.0750 

0.9640 


0.2 190 

0. io c >0 

0. 16 10 

0.0690 

0.0590 


0.2080 

0. 1 790 

0. 1500 

0. 0630 

J . 056' * 


0. 1970 

0.1680 

0. 140 ; 

0.0590 

0.0520 


0. i 860 

0 .13 7 0 

0. 1300 

Q . u 550 

0 .0490 


0. 1 750 

0.147 0 

0. 120 ) 

0.0510 

9.0460 


0.1640 

>. 1 360 

0. 1 100 

0.04/0 

0.0430 


0.1 530 

9. 1240 

0. 1000 

0.0430 

J . 0400 


0.1410 

9.1139 

0. 08 90 

0.0450 

9 . 0 ' / 9 


0. 1 300 

0 . 1 U 2 0 

0. 0790 

-0. 

-9. 


0.1 1 H 0 

0 . GV 1 9 

0.0690 

-0. 

-0. 


0.1 GbO 

0 .0809 

0. 0590 

-0. 

-■J . 


0.09 JO 

) • 0 7 G J 

0. 0500 

- o . 

r> 

> • 


0.0600 

0.6600 

0. 0400 

-0. 

J • 


G . 3000 

0.4000 

0.5000 

1.0000 

1 . 1000 


1.0500 

1 . 1 500 

l .2500 

1.7500 

1 . 8500 

1.9500 
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STAG. SPEED C F SOUND 
ITERATION NO. 1 
1 TEH AT I UN NO. I 
ITERATION NO. 3 
ITERATION NO. 4 
ITERATION NO. 5 
ITERATION NO. 6 
ITERATION NU. 7 
ITERATION NO. 8 
ITERATION NU. 9 
ITERATION NU. 10 
ITERATION NO. LL 
ITERATION NO. IZ 
ITERATION NO. 13 
ITERATION NO. 14 
ITERATION NO. 15 
ITERATION NO. 16 
ITERATION NO. 17 
ITERATION NO. 18 
ITERATION NO. 19 
ITERATION NO. 20 
ITERATION NO. 21 
ITERATION Nu. 22 
ITERATION NO. 23 
ITERATION NO. 24 
ITERATION NO. 25 
ITERATION NU. 26 
ITERATION NO. 27 
ITERATION NU. 28 
ITERATION NO. 29 
ITERATION NU. 30 
ITERATION NU. 31 
ITERATION NO. 32 
ITERATION NO. 33 
ITERATION NU. 34 
ITERATION NO. 35 
ITERATION NO. 3o 
ITERATION NO. 37 
ITERATION NU. 3b 
ITERATION NU. 39 
ITERATION NU. 40 
ITERATION NU. <tl 
ITERATION NO. 42 
ITERATION NO. 43 
I TERAT ION NO . 44 
ITERATION NU. 45 
ITERATION NO. 4o 
ITERATION NU. 47 


AT INLET = 1192.22 

MAX. STREAMLINE 
MAX. S 1 RE AML I NE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. streamline 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
M4A. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 
MAX. STREAMLINE 


CHANGE 

x 

0.2221)3 

CHANGE 

= 

0. 20d lo 3 

CHANGE 

= 

0.164211 

CHANGE 

x 

0. 161841 

CHANGE 


0. 141914 

CHANGE 

= 

0. 124332 

CHANGE 

- 

0. 1 0698 7 

CHANGE 

- 

0.096636 

CHANGE 

- 

0.084035 

CHANGE 

x 

0.0739^2 

CHANGE 

- 

0.06 51 SO 

CHANGE 

= 

0.057468 

CHANGE 

= 

0. ! 50 7 36 

CHANGE 

= 

0.044H2 l 

CHANGE 

x 

0.039616 

CHANGE 

x 

0.035032 

CHANGE 

- 

0. 030968 

CHANGE 

x 

0.027419 

CHANGE 

- 

0.02426 7 

CHANGE 

- 

0.021483 

CHANGE 

x 

0.01902 3 

CHANGE 

= 

0.016847 

CHANGE 

x 

0.014923 

CHANGE 

x 

0.013248 

CHANGE 

= 

0.011771 

CHANGE 

- 

0.0 10461 

CHANGE 

= 

0.009301 

CHANGE 

= 

0.008270 

CHANGE 

x 

0.007356 

CHANGE 

- 

0.006546 

CHANGE 

- 

0.005824 

CHANGE 

- 

0.00521 9 

CHANGE 

- 

0.00.610 

CHANGE 

- 

0.004130 

CHANGE 

x 

0. 0036 2 

CHANGE 

- 

0.003232 

CHANGE 

x 

0.002925 

CHANGE 

= 

0.002615 

CHANGE 

X 

0.002326 

CHANGE 

X 

0.002069 

CHANGE 

= 

0.001845 

CHANGE 

X 

0.001645 

CHANGE 

- 

0. 0)146 7 

CHANGE 

X 

0.09131 5 

CHANGE 

= 

0.001 1 67 

CHANGE 

= 

0.00 104 7 

CHANGE 

- 

0.00094 7 
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T ER AT I ON NO. 48 MAX. STREAMLINE CHANGE “ 0.000833 



ITERATION NO. 49 
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TEKATI0N NO. 49 MAX. STREAMLINE CHANGE = 0.000750 
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